function [xm, fv] = icyPSO(fitness, N, c1, c2, w, t, dim)

  % 1. Determine the number of particle and iteration times
  % 2. Initialize position and velocity of every particle randomly
  % 3. Update position and velocity of every particle using:
  %    V_present = w * V_next + c1 * rand * (P_partOpt - P_next)
  % + c2 * rand * (P_globalOpt - P_next)
  %    P_present = P_next + V_present
  % 
  % fitness     --- {} target function
  %    N        --- {} particle number
  %  c1,c2      --- {} constant
  %    w        --- {} weight of speed
  %    t        --- {} iteration time
  %   dim       --- {} dimension of problem
  %    X        --- { N * dim } position
  %    Y        --- { N * dim } fitness(X)
  %    V        --- { N * dim } velocity
  % X_partOpt   --- { N * dim } part optimization solution of every particle
  % X_globalOpt --- { 1 * dim } global solution of the group
  
  % declare the output number formatter
  format long;

  % initialize position X and velocity V and the other variables
  X = randn(N, dim);
  V = rand(N, dim);

  Y = zeros(N, dim);
  X_partOpt = zeros(N, dim);
  X_globalOpt = zeros(1, dim);

  % initialize fitness results and part optimize result
  for index = 1 : N
    Y(index, :) = fitness(X(index, :));
    X_partOpt(index, :) = X(index, :);
  end
  
  % initialize global optimize result
  X_globalOpt = X(N, :);
  
  for index = 1 : (N - 1)
    if fitness( X(index, :) ) < fitness(X_globalOpt)
      X_globalOpt = X(index, :);
    end
  end

  % interation & update & optimization
  for indexT = 1 : t
    for indexN = 1 : N
      % update velocity
      V(indexN, :) = w * V(indexN, :) + c1 * rand * (X_partOpt(indexN, :) - X(indexN, :)) + c2 * rand * (X_globalOpt - X(indexN, :));

      % update position
      X(indexN, :) = X(indexN, :) + V(indexN, :);

      % optimize
      if fitness( X(indexN, :)) < Y(indexN, :)
	Y(indexN, :) = fitness(X(indexN, :));
	X_partOpt(indexN, :) = X(indexN, :);
      end

      if Y(indexN, :) < fitness(X_globalOpt)
	X_globalOpt = X_partOpt(indexN, :);
      end
    end

    Ybest(indexT) = fitness(X_globalOpt);
  end
  
  xm = X_globalOpt';
  fv = fitness(X_globalOpt);
end

